################################################################################
### Technological Change and the International System Replication Files
### Figure 3. Substantive Effects - The International System and Technology Adoption
###  
### Required data files: "tech_data.dta"
###                    
### Created by: Thomas Cunningham
### Date: 27 January 2021
##################################################################

#Start from Clean Work Space
rm(list = ls())

#Set Working Directory
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) 
# Note: if you are not using R Studio this command will not work, set WD to source file location manually
require(readstata13)
data <- read.dta13("tech_data.dta")


delta <- TRUE
library(Zelig)


for(iv in c("sysconall", "nviable", "c4")[1]){
  
  if(delta == TRUE){
    temp_dat <- na.omit(data[, c("delta", 
                                 "log_adopt", 
                                 iv, 
                                 "polity2", 
                                 "warL1", 
                                 "cwarL1", 
                                 "linear_time_trend", 
                                 "cou_tech_enc", 
                                 "sdt")])
    
    colnames(temp_dat)[3] <- "my.iv"
    
    model <- zelig(delta ~ my.iv + sdt + polity2 + warL1 + cwarL1, 
                   data = temp_dat,
                   model = "normal",
                   cite = FALSE)
    
    # Set number of simulations
    sim.n <- 5000
    
  } else {
    print(Sys.time())
    
    print("Fitting model")
    
    temp_dat <- na.omit(data[, c("log_adopt", 
                                 iv, 
                                 "polity2", 
                                 "warL1", 
                                 "cwarL1", 
                                 "linear_time_trend", 
                                 "cou_tech_enc", 
                                 "sdt")])
    
    colnames(temp_dat)[2] <- "my.iv"
    
    
    model <- zelig(log_adopt ~ my.iv + sdt + polity2 + warL1 + cwarL1 + linear_time_trend + as.factor(cou_tech_enc), 
                   data = temp_dat, 
                   model = "normal", 
                   cite = FALSE)
    
    print("Model fit completed.")
    
    sim.n <- 10
  }
  
  model.variables <- c("my.iv", 
                       "sdt", 
                       "polity2", 
                       "warL1", 
                       "cwarL1")
  
  iv.var.name <- c("Syscon (Singer 1972)", 
                   "Viable Great Power Coalitions",
                   "C4")[which(c("sysconall", 
                                 "nviable", 
                                 "c4") == iv)]
  
  variable.names <-  c(iv.var.name, 
                       c("Spatial Distance to \nTechnology, 3 year lag", 
                         "Polity2 score", 
                         "War (lagged one year)", 
                         "Civil War (lagged one year)", 
                         "Linear Time Trend"))[1:length(model.variables)]
  
  name.ind <- 1
  
  for(variable in model.variables){
    range <- rep(mean(temp_dat[, variable]), 3)
    range[1] <- range[1] - sd(temp_dat[, variable])
    range[3] <- range[3] + sd(temp_dat[, variable])
    # }
    
    pdata <- data.frame(matrix(ncol = length(range), nrow = sim.n))
    colnames(pdata) <- c("down", "mean", "up")
    
    pdata$var <- variable.names[name.ind]
    name.ind <- name.ind + 1
    
    print(paste("Simulating -1SD, MEAN, +1SD for", variable.names[name.ind-1]))
    
    index <- 0
    for(i in 1:length(range)){
      
      # Set level of IV  
      model.fit <- eval(parse(text = paste("setx(model, ", variable, " = ", range[i], ")")))
      
      # Simulate quantities of interest
      sim   <- sim(model, x = model.fit, num = sim.n)
      
      # Get and save expected values
      pdata[, i] <- unlist(sim$get_qi(qi = "ev", xvalue = "x"))
      
      cat(paste("..", i))
      
    }
    
    cat("\n")
    
    big.pdata <- if(exists("big.pdata")){
      rbind(big.pdata, pdata)
    } else {
      pdata
    }
    
  }
  
  print("Simulation complete.")
  
}

big.pdata$var <- factor(big.pdata$var, 
                        levels = rev(c("Syscon (Singer 1972)", 
                                       "Spatial Distance to \nTechnology, 3 year lag", 
                                       "Polity2 score", "War (lagged one year)", 
                                       "Civil War (lagged one year)", 
                                       "Linear Time Trend")))


library(ggplot2)


bpdata <- subset(big.pdata, 
                 big.pdata$var != "Linear Time Trend")

bpdata_full <- bpdata

bpdata <- bpdata_full
bpdata$mid <- mean(bpdata$mean)
bpdata <- bpdata[!bpdata$var == "Spatial Distance to \nTechnology, 3 year lag", ]
bpdata$down_mean <- ave(bpdata$down, bpdata$var, FUN = function(x) mean(x))
bpdata$up_mean <- ave(bpdata$up, bpdata$var, FUN = function(x) mean(x))
bpdata$down_95 <- ave(bpdata$down, bpdata$var, FUN = function(x) {sort(x)[ceiling(length(x)*0.95)]})
bpdata$down_05 <- ave(bpdata$down, bpdata$var, FUN = function(x) {sort(x)[floor(length(x)*0.05)]})
bpdata$up_95 <- ave(bpdata$up, bpdata$var, FUN = function(x) {sort(x)[ceiling(length(x)*0.95)]})
bpdata$up_05 <- ave(bpdata$up, bpdata$var, FUN = function(x) {sort(x)[floor(length(x)*0.05)]})
bpdata$jitterdown <- bpdata$down
bpdata$jitterdown[bpdata$jitterdown >= bpdata$down_95] <- NA
bpdata$jitterdown[bpdata$jitterdown <= bpdata$down_05] <- NA
bpdata$jitterup <- bpdata$up
bpdata$jitterup[bpdata$jitterup >= bpdata$up_95] <- NA
bpdata$jitterup[bpdata$jitterup <= bpdata$up_05] <- NA


for(j in unique(bpdata$var)){
  if(mean(bpdata$jitterdown[bpdata$var == j], na.rm = TRUE) < mean(bpdata$jitterup[bpdata$var == j], na.rm = TRUE)){
    bpdata$down[bpdata$var == j] <- NA
    bpdata$down_mean[bpdata$var == j] <- NA
    bpdata$jitterdown[bpdata$var == j] <- NA
    bpdata$down_05[bpdata$var == j] <- NA
    bpdata$down_95[bpdata$var == j] <- NA
    
  } else {
    bpdata$up[bpdata$var == j] <- NA
    bpdata$up_mean[bpdata$var == j] <- NA
    bpdata$jitterup[bpdata$var == j] <- NA
    bpdata$up_05[bpdata$var == j] <- NA
    bpdata$up_95[bpdata$var == j] <- NA
  }
}

p <- ggplot(unique(bpdata), aes(x=mid*100, y = var, group = var))+
  geom_vline(xintercept=bpdata$mid[1]*100, linetype = 3)+
  geom_jitter(aes(x=jitterdown*100), height = 0.05, alpha = 0.05, col = "gray19")+
  geom_jitter(aes(x=jitterup*100), height = 0.05, alpha = 0.05, col = "gray65")+
  geom_point(aes(x=down_mean*100), shape = 16, size = 4, colour = "white")+
  geom_point(aes(x=up_mean*100), shape = 17, size = 4, colour = "white")+
  geom_point(aes(x=down_mean*100, col = "-1 SD"), size = 3, shape = 16, alpha = 1)+
  geom_point(aes(x=up_mean*100, col = "+1 SD"), shape = 17, size = 3)+
  theme_bw()+
  xlab("% Change in Technology Units per Capita), per year")+
  ylab("")+
  theme(legend.title = element_blank(), legend.pos = "bottom", plot.margin = margin(0, 0, 0, 0, "cm"))+
  scale_colour_manual(breaks=c("-1 SD","+1 SD"), name='',values=c("-1 SD"="gray19", "+1 SD"="gray65"))+
  guides(colour=guide_legend(override.aes=list(shape=c(16, 17), linetype = c("blank", "blank"))))+
  scale_shape_manual(values=c(16, 17))
p

p <- p+scale_x_continuous(breaks = c(bpdata$mid[1]*100,4.5,5,5.5,6), labels = c("4.26\n(Mean)","4.5","5","5.5","6"))

ggsave(file = "output/figure3.pdf", width=6, height=3, device = "pdf")
